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Hydrodynamics  of  Gas  Channel  Formation 


1.  INTRODUCTION 

There  has  been  interest  recently  in  the  formation  of  heated,  reduced-density  channels  in 
gaseous  media,  by  passing  a  laser  pulse,1  an  electron  beampulse,2  an  electric  discharge,3  or  an 
exploding  wire  current4  through  the  gas.  Such  a  channel  may  be  used  to  guide  a  charged  particle 
beam  from  the  beam  source  to  a  fusion  pellet  in  an  inertial  confinement  fusion  reactor.  In 
some  cases,  it  may  also  be  possible  to  pulse  the  channel  several  times  by  passing  widely  spaced 
energy  pulses  through  it. 

In  this  report,  we  study  specifically  the  hydrodynamics  of  channel  formation.  Analytic 
solutions  have  been  obtained  in  the  past  for  blast  wave  hydrodynamics  in  cylindrical  geometry,5 
but  only  for  strong  blast  waves  (where  the  shock  has  high  Mach  number,  even  after  it  has  ex¬ 
panded  to  many  times  the  initially  heated  radius).  Such  solutions  are  valid  only  if  the  initial 
overpressure  is  hundreds  or  thousands  of  times  the  ambient  pressure.  Our  interest  here  is  in 
more  moderate  overpressures.  One-dimensional  fluid  code  simulations  of  the  radial  gas  hydro¬ 
dynamics  are  presented,  over  a  wide  range  of  parameters.  Theoretical  arguments,  as  well  as 
scaling  laws  induced  from  the  simulations,  are  used  to  construct  analytic  models  which  include 
all  the  qualitative  features  of  channel  hydrodynamics,  and  permit  accurate  quantitative  calcula¬ 
tion  of  the  process  without  further  resort  to  hydrodynamic  codes.  In  order  to  focus  as  clearly  as 
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possible  on  the  hydrodynamics,  we  avoid  getting  into  beam  propagation  and  energy  deposition, 
late-time  channel  cooling,  and  the  complications  of  real  gas  thermal  properties. 

We  consider  the  case  in  which  the  duration  of  the  heating  pulse  is  so  short  that  it  can  be 
treated  as  instantaneous  on  the  hydrodynamic  time  scale.  In  the  case  of  multiple  heating 
pulses,  we  treat  the  regime  in  which  pulse  spacing  is  long  enough  so  that  low  density  gas  chan¬ 
nel  will  have  (or  nearly  have)  come  to  rest,  in  pressure  balance  with  the  ambient  gas,  by  the 
time  the  next  pulse  arrives.  Thus  the  hydrodynamic  process  is  an  interative  one,  in  which  each 
pulse  serves  only  to  leave  behind  an  essentially  static  channel,  i.e.,  an  initial  condition  for  the 
next  stage  of  channel  expansion.  Thus  the  problem  we  set  for  ourselves  is  as  follows.  Given  a 
gas  density  profile  /„  =  0),  created  by  the  first  n  —  1  pulses;  t„  is  time  measured  with 

(«'0  at  the  arrival  time  of  the  nlh  pulse.  Given  also  a  pressure  profile  P„(r,  t„  =  0+)  = 
P„-\(r)  +  P„(r),  where  /Vilr),  which  is  close  to  the  ambient  pressure,  is  due  to  the  previous 
pulses,  and  Pn(r)  is  the  overpressure  profile  due  to  energy  deposition  at  tn  —  0  by  the  n'h  pulse. 
We  calculate  the  subsequent  hydrodynamic  response  of  the  gas. 

The  properties  of  the  heating  source  (laser,  electron  beam,  or  whatever)  appear  in  our  cal¬ 
culations  only  through  the  "initial"  overpressure  profile,  after  the  passage  of  the  n  th  pulse, 

P„(r).  For  specificity  we  assume  a  Bennett  profile,  P„(r)  =  - - -  ,  a  typical  smooth 

(1  + 

bell-shaped  profile  which  is  particularly  appropriate  for  an  electron  beam  under  certain  condi¬ 
tions.6  We  shall  extend  this  work  to  consider  arbitrary  profiles  in  a  subsequent  paper. 

A  further  simplification  is  made  in  the  presentation  here,  by  taking  the  adiabatic  index  of 
the  gas,  y  =  Cp/ Cv,  to  be  constant  (where  Cp  and  Cv  are  the  specific  heats  at  constant  pressure 
and  volume).  In  numerical  work,  we  assume  specifically  y  =  7/5,  the  value  appropriate  to  a 
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diatomic  gas  without  internal  degrees  of  freedom.  In  real  gases,  y  can  vary  considerably  at  high 
temperatures,  where  internal  degrees  of  freedom  are  excited;  such  complications  will  be  includ¬ 
ed  in  future  work. 

The  principal  results  of  this  report  will  be  summarized  in  the  remainder  of  this  section. 
Section  2  presents  the  detailed  simulation  and  theory  for  a  first  pulse,  fired  into  a  homogeneous 
gas.  Section  3  gives  analogous  results  for  a  pulse  depositing  its  energy  in  a  preformed  channel. 
Section  4  deals  with  the  special  case  of  a  weak  pulse,  which  results  in  an  overpressure  small 
compared  to  the  ambient  pressure;  in  this  case,  the  channel  formation  process  reduces  to  linear 
sound  wave  propagation,  and  the  complete  time-dependent  hydrodynamics  can  be  solved  in 
closed  form.  Section  5  briefly  summarizes  our  conclusions. 

In  order  to  elucidate  the  scaling  laws,  it  is  best  to  deal  with  the  channel  formation  process 
in  dimensionless  variables,  scaled  to  ambient  values,  with  r  scaled  to  heating-source  (or 
"beam")  radius  a ,  and  t  scaled  to  a/cs,  wherer  cs  is  ambient  sound  speed.  These  dimensionless 
variables  are  denoted  by  tildes.  Table  1  summarizes  the  notation  used  throughout  the  paper. 

The  hydrodynamic  processes  following  the  passage  of  a  beam  pulse  with  a  Bennett  profile 
can  be  separated  into  three  phases,  on  different  time  scales.  (1)  First,  gas  density  reduction  by  a 
rarefaction  wave  within  the  beam  radius  (2)  Expansion  and  density  reduction  within  the 
rarefied  channel,  which  is  typically  several  times  larger  than  the  beam  radius.  In  general,  the 
channel  as  a  whole  reaches  a  minimum  density  state,  when  the  pressure  throughout  it  is  (over  a 
wide  range  of  parameters)  about  85-90%  of  ambient  pressure.  (3)  Slow  return  of  the  channel 
to  ambient  pressure,  with  an  accompanying  density  increase  of  about  10%.  The  time  scales  for 
these  processes  are  summarized  in  Table  2. 
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If  the  pulse  separation  is  long  enough  so  that  the  channel  reaches  a  stationary  state 
between  pulses,  the  hydrodynamic  process  is  an  iterative  one,  in  which  the  density  profile  p„(r ) 
after  n  pulses  can  be  expressed  in  terms  of  the  previous  channel  profile,  p„_,(r),  and  the  over¬ 
pressure  on  axis  due  to  the  n'h  pulse,  Pon.  The  density  reduction  at  r  —  0  (and,  in  general,  for 
r  <  a)  is  due  simply  to  adiabatic  expansion,  and  can  be  expressed  in  a  very  simple  analytic 
form.  In  fact,  adiabatic  expansion  is  a  fairly  reasonable  model  for  the  channel  profile  as  a 
whole,  but  the  density  channel  is  widened  by  up  to  -20%  per  pulse  by  additional  heating,  due 
to  a  radially  outward-propagating  shock  driven  by  each  heating  pulse.  The  results  for  the  den¬ 
sity  profile  reached,  when  the  channel  comes  to  rest  after  the  n'h  pulse,  are  summarized  in 
Table  3. 

If  the  pulses  are  spaced  such  that  each  one  arrives  when  the  channel  is  at  minimum  densi¬ 
ty  (at  pressure  10-15%  below  ambient),  then  the  density  profiles,  at  the  time  when  the  n'h  pulse 
arrives,  are  simply  depressed  below  the  results  of  Table  3  by  about  10%.  This  effect  does  not 
increase  cumulatively  from  pulse  to  pulse. 

Finally,  we  note  that  for  an  initial  Bennett  overpressure  profile,  the  channel  does  not  un¬ 
dergo  repeated  strong  radial  oscillations  after  heating  by  any  given  pulse.  Rather,  its  density 
dips,  only  once,  about  10%  below  its  final  value,  and  then  slowly  returns  to  its  final  value.  If 
"ringing"  of  the  channel  occurs,  it  is  at  very  small  amplitude. 

Figures  1.1  through  1.3  give  a  quick  sample  of  some  results  derived  by  the  analytic  argu¬ 
ments  described  above.  Fig.  1.1  shows  the  central  density  when  the  gas  comes  to  rest  after  the 
first  pulse,  as  a  function  of  the  central  overpressure  P0.  Figure  1.2  shows  the  coordinate  r,  of 
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the  fluid  element  that  originated  at  r  =  1,  when  the  gas  comes  to  rest  after  the  first  pulse.  Fig¬ 
ure  1.3  shows  a  graphical  method  of  calculating  the  location  of  any  given  fluid  element,  when 
the  gas  comes  to  rest  after  the  first  and  second  pulses. 

II.  GAS  EXPANSION  FOLLOWING  PASSAGE 

OF  A  SINGLE  PULSE 

A.  General  Remarks 


In  this  section  we  present  extensive  numerical  solutions  for  the  radial  expansion  of  the 
gas  following  the  passage  of  a  single  pulse,  outline  the  principal  features  of  this  hydrodynamic 
evolution,  and  present  analytical  and  heuristic  techniques  that  treat  all  of  these  features  and 
permit  quantitative  calculation  of  the  evolution  by  much  simpler  methods.  For  simplicity,  we 
assume  a  constant  value  for  the  ratio  of  specific  heats  y  =  cp/cv ,  and  choose  y  =  7/5,  as  is  ap¬ 
propriate  for  a  diatomic  gas  in  the  absence  of  internal  degrees  of  freedom.  We  assume  that  the 
gas  is  initially  heated  by  the  (essentially  instantaneous)  passage  of  the  pulse,  and  that  the  initial 
heating  profile  is  the  Bennett  profile,  i.e.,  the  initial  conditions  for  density  p,  pressure  P,  and 
temperature  T  are 


P(r,  0+) 


p(r,0)  =  pA. 


- Pa 


(l+;'/.a?)Jr 


(2.1) 

(2.2) 


T(r,  0)  -  Ta  ryy  (2.3) 

where  the  subscript  A  denotes  ambient  values.  Thermal  conduction  and  radiative  cooling  are 
neglected,  since  they  proceed  much  slower  than  the  hydrodynamic  evaluation.  The  only  dissi¬ 
pation  included  in  the  model  is  that  which  occurs  in  shocks;  the  latter  is  treated  in  a  manner 
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which  reproduces  the  Rankine-Hugoniot  relations  with  great  accuracy,  without  resorting  to  the 
(sometimes  troublesome)  numerical  artifice  of  a  Von  Neumann  viscosity. 

It  is  natural  to  treat  this  problem  in  dimensionless  units 

~r  =  ''/a,  (2.4a) 

<  “  tcA/a,  (2.4b) 

P=  Pip  a-  P  =  PI  PA  (2.5) 

where  a  is  the  Bennett  radius  and  cA  =  (yPA/pA)in  is  the  ambient  sound  speed.  In  terms  of 

these  units,  the  evolution  depends  only  on  the  single  dimensionless  parameter  P0>  the  initial 

overpressure  at  r  =  0,  scaled  to  the  ambient  pressure.  For  example,  if  the  Bennett  radius  a  is 

halved,  but  the  total  energy  deposited  is  held  constant,  then  P0  is  quadrupled. 

This  initial  value  problem  has  been  solved  numerically  for  a  wide  range  of  values  of  P0. 
The  one-dimensional  (radial)  code  ETBFCT  was  used  to  solve  the  Euler  hydrodynamic  equa¬ 
tions  (conservation  equations  of  mass,  momentum  and  energy)  in  conservative  form.  The  code 
is  Eulerian,  with  non-uniform  grid  spacing  and  regridding  capability  to  enhance  accuracy  near 
the  origin  and  at  the  shock  location.  The  code  uses  the  flux-conserving  transport  (FCT)7  tech¬ 
nique,  and  has  been  tested  extensively  to  verify  its  accuracy  in  treating  shocks  in  non-uniform 
media."  It  has  been  shown  to  reproduce  the  jump  conditions  across  a  discontinuity  (i.e.  the 
Rankine-Hugoniot  conditions  across  a  shock)  with  extraordinary  accuracy,  if  one  looks  several 
cells  downstream  from  the  shock.  However,  some  averaging  occurs  over  two  or  three  cells  at 
the  shock,  which  numerically  clips  the  density  and  pressure  peak  at  the  shock,  preventing  one 
from  accurately  reading  off  the  shock  strength  in  a  strongly  non-uniform  case.  This  effect  is 
reduced  (to  a  factor.<1.2  in  our  computer  runs)  by  using  a  fine  grid  spacing  near  the  shock.  In 
the  various  figures  showing  density  and  pressure  profiles,  we  have  extrapolated  the  numerical 
results  to  reconstruct  sharp  peaks  at  the  shocks. 
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The  time  evolution  of  the  density  and  pressure  profiles  is  shown  in  Figs.  2.1-8,  for  cases 
with  P0  =  42.0,  10.5,  1.8,  and  0.1.  In  all  cases,  the  following  sequence  of  events  occurs,  (i)  A 
pressure-density  pulse  forms  on  the  outside  of  the  heated  gas  region,  and  begins  to  move  out 
radially.  A  rarefaction  forms  behind  this  pressure  pulse,  causing  the  central  density  to  drop, 
(ii)  For  a  smooth  (e.g.  Bennett)  initial  pressure  profile,  the  density  and  pressure  pulses  are  ini¬ 
tially  smooth,  but  steepen  to  a  shock  wave  after  a  finite  time*  which  will  be  estimated  theoreti¬ 
cally  in  Sec.  2C.  1  The  shock  strength  (indicated,  e.g.,  by  the  density  jump)  grows  for  a  while, 
reaches  a  maximum,  and  then  faiis  off  because  of  the  cylindrical  expansion,  (iii)  Behind  the 
shock,  a  low-density  channel  forms,  the  temperature  decreases  somewhat  from  its  initially  high 
value  because  of  PdV  work  done  against  the  ambient  presure,  and  the  pressure  falls  until  it 
reaches  the  ambient  level,  (iv)  The  shock  eventually  detaches  from  the  channel,  and  propagates 
toward  r  =  °°,  carrying  part  of  the  initially  deposited  energy  with  it.  (v)  Meanwhile,  the  density 
in  the  channel  undershoots  pressure  equilibrium  with  the  ambient  gas  by  a  small  amount,  and 
then  the  density  and  pressure  increase  very  slowly  until  pressure  equilibrium  is  reached.  The 
density  channel  does  not  undergo  a  series  of  violent  oscillations  about  pressure  equilibrium  such 
as  occurs  in  the  well-known  case  of  an  underwater  explosion,9  i.e.,  a  high-pressure  gas  bubble 
formed  in  an  essentially  incompressible  medium.  On  the  contrary,  only  a  single,  weak 
overshoot  of  equilibrium  occurs. 

B.  Final  State 

The  density  profile  at  late  times,  after  the  gas  has  come  to  rest  in  pressure  equilibrium, 
can  be  expressed  in  terms  of  the  initial  state  and  the  properties  of  the  shock.  The  specific  en¬ 
tropy  of  a  gas  is  expressed  in  terms  of  y  as 


*ln  the  cases  of  weak  initial  heating,  the  shock  is  so  weak  that  it  is  difficult  to  resolve,  and  is  of  no  real  importance. 
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5  =  cv  In  ( Pp  y).  (2.6) 

In  our  system,  dissipation  occurs  only  in  the  shock  front.  Any  given  element  of  fluid  is  over¬ 
run  by  the  shock  no  more  than  once.  If  we  let  A S(r0)  be  the  specific  entropy  increase  when 
the  shock  crosses  the  fluid  element  whose  initial  position  was  r0,  then  Eq.  (2.6)  can  be  rewritten 


P(r  ,t)  Il/y 

~p{r°'l)~  Tb°, 0+)~  exp[-AS(r0)/cpl,  (2.7) 

where  we  use  the  dimensionless  variables  (2.5).  In  Eq.  (2.7),  is  used  as  a  Lagrangian  coordi¬ 
nate,  i.e.,  a  label  for  the  fluid  element  located  at  r0  at  t  -  0.  The  Eulerian  coordinate  r,  i.e„ 
the  position  of  this  fluid  element  at  time  /,  can  be  regarded  as  a  function  of  r0  and  /.  At  t  =  °o 


in  particular, 


Ptf,.™)-  1, 


p(?0,°o)  =  [P(?o,0+l]-^  exp  [-A5(?0)/c.] 


~  I1  +  Tr+r^J  exp  hAS(?o)/cp]-  (2.8) 

In  (2.8),  we  use  also  a  dimensionless  Lagrangian  coordinate, 

r0~rja.  (2  9) 

Since  Eqs.  (2.7),  (2.8)  are  written  in  terms  of  the  Lagrangian  coordinate  r0,  the  Eulerian  coor¬ 
dinate  r(r0,t )  must  be  prescribed  to  complete  the  description.  This  is  accomplished  by  noting 
that  the  volume  element  d(irr2)  is  inversely  proportional  to  the  density  p,  so  that 


-ii-  f2{r  t)  -  ?(?o,0)  -  1 

ft!  p(r0,~t)  p  (?0j ) 


(2.10a) 


f2(f0,t)  -  d(r0)2/p(?0,t). 


(2.10b) 
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or,  in  particular,  at  late  times, 

r2(r0,°°)  =  JQ  °  d(r0)2lp(r0,<x>).  (2.10c) 

Equations  (2.8)  and  (2.10c)  completely  specify  the  density  profile  at  t  =  °°,  once  AS  is  known. 


Since  the  shocks  are  rather  weak,  even  for  overpressures  up  to  P  =  42  (as  will  be  dis¬ 
cussed  in  Sec.  2c),  the  approximation 


AS  =  0  (2.11) 

is  qualitatively  reasonable.  (We  shall  see  that  it  is  exact  within  the  central  channel.)  If  this  ap¬ 
proximation  is  made,  Eq.  (2.8)  reduces  to 


p(,r0,°°)  =  [P(ro,0+)l  1/y  = 


1  + 


(1  +  r2)2 


-l/y 


(2.12) 


and  Eqs.  (2.10c)  and  (2.12)  give  a  very  simple  first  approximation  to  p(r,°°).  This  result  is 
plotted  as  solid  curves  in  Figs.  2.9-2.12,  for  the  various  values  of  P0,  along  with  the  computed 
density  profiles  at  late  time.  We  note  that  the  approximation  of  Eqs.  (2. 10c, 12)  predicts  the 
depth  of  the  central  density  hole  exactly,  after  the  gas  has  come  to  pressure  equilibrium.  How¬ 
ever,  the  approximation  underestimates  the  width  of  the  density  channel  by  up  to  20%  in  the 
strongest  case.  Of  lesser  importance,  the  approximation  fails  to  predict  a  broad,  very  shallow 
(<  5%)  wing  to  the  density  channel.  As  will  be  discussed,  the  channel  is  wider,  but  not 
deeper,  than  the  isentropic  model  prediction,  because  only  the  outer  regions  are  shock-heated. 


C.  Shock  Formation  and  Entropy  Increase 

We  note  first  that  the  radially  outward-propagating  shock  forms  at  some  radial  position  rs , 
and  overruns  all  of  the  gas  outside  rs.  The  central  region  r  <  rs  is  never  overrun  by  a  shock, 
and  its  expansion  is  thus  ijentropic,  and  given  by  Eqs.  (2.10c)  and  (2.12),  in  agreement  with 
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the  fluid  code  results.  The  code  runs  indicate  that  the  shock  forms  in  the  fluid  element  whose 
initial  (i.e.,  Lagrangian)  coordinate  is 

h  =  hs  —  0.8;  (2.13) 

this  will  shortly  be  explained  theoretically.  Thus  the  entire  central  region  of  the  low  density 

channel  (e.g.,  out  to  r  —  3  after  expansion  has  occurred,  for  P0  =  42)  should  be,  and  is,  accu¬ 
rately  predicted  at  t  =  °°  by  Eqs.  (2.10c)  and  (2.12). 


We  next  discuss  the  determination  of  r05,  in  a  qualitative  way.  Any  hydrodynamic  distur¬ 
bance  will  eventually  steepen  and  form  a  shock,  but  the  time  ts  necessary  for  this  to  occur 
varies  inversely  with  the  strength  of  the  disturbance.  We  know  of  no  general  formula  for  is 
(except  in  the  idealized  case  of  a  Riemann  simple  wave)  but  we  can  construct  an  approximate 
formula  as  follows.  An  idealized  pressureless  fluid  free  streams  according  to  its  initial  velocity 
distribution,  i.e  ,  satisfies  the  equation 


or  equivalently, 


du  +  dw  =  OP 

dt  dr  p  dr 


(2.14a) 


u(r,  t)  =  u(r  -  u(r,  t)t). 


(2.14b) 


It  can  be  shown  exactly10  that,  according  to  this  model,  a  shocklike  discontinuity  forms  in  each 

fluid  element  whose  initial  location  ros  is  such  that  is  a  mjnjmum  an(j  at  tjme 

dr„ 


du(rn  0) 

l - ~~ 

dr0 

However,  in  the  real  problem  under  consideration. 


(2.15) 

the  initial  velocity  is  zero  everywhere,  and  it 


takes  a  finite  time  for  the  velocity  field  to  develop;  also,  the  velocities  continue  to  change  under 
the  influence  of  pressure  gradients.  But  if  the  shock  forms  before  the  pressure  profile  is  greatly 
changed,  which  occurs  if  the  pressure  perturbation  is  strong,  then  it  is  reasonable  qualitatively 
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to  picture  shock  formation  as  a  two-stage  process:  during  a  period  r(,  a  velocity  profile  «(/■) 
develops  out  of  the  initial  pressure  profile  P(r,  0),  and  during  a  subsequent  time  r2,  this  veloci¬ 
ty  profile  free-streams  until  it  steepens  to  a  shock.  The  shock  formation  time  is  then 

ts  —  T)  +  t2.  (2.16) 

At  time  r,,  according  to  this  picture, 


__  fT i  dt  bP(r,  t)  __  t,  8/>(r,0+) 


,  \  __  f  1  dt  aP(r,  t)  __ 

u(r,  r,)  - - I - r - 

p  Qr 

and  according  to  Eq.  (2.15),  r2  is  given  by 


P  c tr  8r 


__  ...  ,  dwCr.ri)  __  p^  ,82/>(r,0+) 
t2  ~  Min,  | - ^ - 1  ~  Min,  —  | - ~2 - 1 


(2.17) 


(2.18) 


where  the  notation  min,  devotes  the  minimum  with  respect  to  r.  But  r,  is  arbitrary  in  this  pic¬ 
ture;  therefore  one  may  choose  tj  such  as  to  minimize  ts,  i.e.. 


...  ...  I  ,  Pa  i  d2P(r,0+)  j_jl 

ts  -  MinT]  Min, jr,  +  — | - ^ - 1  j 


Performing  the  minimization  over  rt,  we  find 

r  ,-pY2 


dr 


-1/2 


and 


ts  -  2 p|2 


Max, 


d2P(r,  0+) 
dr2 


-1/2 


The  shock  forms  in  the  fluid  element  whose  initial  position  ros  is  determined  by 

d3P(ro,0) 


dr3 


I.  -  o. 


For  the  assumed  Bennett  profile  of  P(r0, 0),  we  find 


(2.19) 


(2.20a) 


(2.20b) 


(2.21) 


r0J-  (3/5) 1/2  «  0.78, 
l  -  2.14  P~U2. 


(2.22a) 

(2.22b) 


LAMPE,  SZU,  AND  KAINER 


These  estimates  of  ts  and  ros  are  found  to  be  in  reasonable  agreement  with  the  fluid  code  results 
P0  >  5,  as  seen  in  Figs.  I  and  2.  For  the  weaker  overpressures,  the  actual  shock  onset 
time  can  be  much  later  than  these  estimates,  but  in  these  cases  the  shock  is  so  weak  as  to  be 
insignificant  anyway. 


When  it  first  forms  at  time  rs,  the  shock  is  very  weak.  As  a  result,  there  is  no  discernible 
discontinuity  in  the  final  density  profile  p(r0,  °o)  between  the  unshocked  region  ra  <  ros  and 
the  shocked  region  r0  >  ros.  After  its  initial  formation,  the  shock  propagates  outward  and  con¬ 
tinues  to  grow,  as  additional  fluid  piles  up  behind  it.  It  reaches  a  maximum  amplitude  at 
Lagrangian  coordinate  rom,  and  then  its  amplitude  decreases,  as  its  area  increases  because  of 
cylindrical  expansion.  Theoretical  results  for  the  variation  of  shock  amplitude  exist  in  the 
literature  for  the  self-similar  (blast  wave)  expansions  that  occur  for  very  strong  overpressures,5 
but  these  results  apply  for  Mach  numbers  M  »  1,  for  which  the  density  jump  p"/p+  across 
the  shock  approaches  six.  They  are  not  particularly  valid  even  for  the  strongest  overpressure 
studied  here,  P  =  42,  for  which  M  <  2.4  and  p~/p+  <  3.3.  As  discussed  in  Sec.  IIA,  it  is  also 
difficult  to  plot  the  shock  amplitude  evolution  accurately  from  the  fluid  code  results.  However, 
a  heuristic  estimate  of  shock  strength  has  been  found,  which  agrees  with  the  code  results  to 
within  their  accuracy,  for  cases  where  the  shock  is  strong  enough  to  be  significant  (P0  >  5). 
This  permits  accurate  calculation  of  the  final  density  profile,  and  correctly  represents  the  shock 
onset  from  Eq.  (2.22a),  and  the  fact  that,  for  r0  large,  shock  stength  should  depend  on  P0  and 
r0  only  through  the  combination  P„r0  2,  i.e.,  the  ratio  of  overpressure  energy  P„a2  to  ambient 
thermal  energy  PAr2.  It  is  convenient  to  express  the  shock  strength  in  terms  of  the  density 
jump  p_/p+  across  the  shock: 


p  /p+  = 


1  r0  <  0.8 

1  +  0.33  P0V 4  f0,  0.8  <  r0  <  3 

1  +  1.7  PV*  r-''2,  3  <  r„ 


(2.23) 
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The  entropy  jump  A s(r0)  and  the  Mach  number  M  are  related  to  p7p+  through 

|V7 


exp(-AS/c.) 


M 2 


6  ~  P  Ip*  1 

r\  o* 

6p7p+  -  1 1 

II 

■o’ 

1 

+ 

Os 

1 

•O 

1 

+ 

(2.24) 

(2.25) 


The  entropy  jump  A S(r0)  from  Eqs.  (2.23),  (2.24)  can  be  used  with  Eqs.  (2.8),  (2.10c)  to  pro¬ 
vide  a  complete  calculational  model  for  the  final  density  profile,  when  the  channel  reaches  pres¬ 
sure  equilibrium.  It  can  be  seen  in  Figs.  2.9  and  2.10  that  the  agreement  with  the  fluid  code 
results  is  very  good. 


D.  Time  Evolution 


Thus  far,  we  have  been  concerned  primarily  with  obtaining  accurate  estimates  of  the  den¬ 
sity  profile  p(r,  °°)  at  late  times,  when  the  pressure  P(r,  <x>)  has  returned  to  ambient.  In  this 
section,  we  shall  consider  the  time  scales  for  radial  expansion,  and  the  minimum  density  that 
occurs  during  the  evolution.  It  is  difficult  to  obtain  theoretical  results  in  this  area,  but  analysis 
of  the  fluid  code  solutions  yields  some  simple  heuristic  generalizations. 


We  note,  from  Figs.  2. 1-2.4,  that  the  density  in  the  channel  reaches  a  minimum  during 
the  evolution,  and  then  slowly  returns  to  its  final  equilibrium  value.  When  viewed  in  terms  of 
the  pressure  profile,  as  in  Figs.  2. 5-2. 8,  this  undershoot  takes  a  very  simple  form  in  the  whole 
range  P0  >  1.  In  all  of  these  cases,  the  pressure  profile  P(r,  t)  takes  on  a  very  broad,  shallow 
minimum  value  Pm,  slightly  below  PA,  at  some  time  tm.  Since  P(r)  is  practically  constant  over 
the  density  channel,  the  density  profile  at  tm  is  given  by  Eqs.  (2.8),  (2.10c),  simply  by  replacing 
PA  with  P„,  i.e.. 


P(r„.  O 


51/2 


1  + 


(i  +  h2)2 


-l/y 


exp  [-  ASCr0)/c/>], 


(2.26) 
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and 

'2('o.  U  -  f0  °  d(r0)2/  p(r0,  tj.  (2.27) 

The  pressure  minimum  Pm  is  plotted  as  a  function  of  P0  in  Fig.  13.  It  is  seen  to  be  very  slowly 
varying  over  the  whole  range  P0  >  1;  one  is  not  far  off  in  estimating  P„  -  0.88  for  all  such 
values  of  P0.  Thus,  in  Lagrangian  coordinates,  the  minimum  density  p(r0,  tm)  is  approximately 
10%  less  than  the  equilibrium  profile  p(T0,  °°).  In  Eulerian  coordinates,  the  effect  is  greater, 
since  the  channel  itself  is  also  widened  by  approximately  a  factor  (1.10) 1/2  ~  1.05. 

This  particularly  simple  picture  cannot  hold  in  the  regime  P0«  1,  where  all  perturbed 
quantities,  e.g.,  1  -  Pm,  must  be  linearly  proportional  to  P0  (as  shown  in  Fig.  2.13).  If 
P 0  «  1  the  channel  formation  process  becomes  one  of  linear  sound  wave  propagation.  In 
such  a  process,  the  pressure  profile  is  of  similar  width  (rather  than  much  broader)  to  the  densi¬ 
ty  profile.  A  complete  time-dependent  analytic  solution  for  the  linear  case  P  «  1  is  given  in 
Sec.  IV. 

The  time  evolution  of  the  density  channel  can  be  divided  roughly  into  three  stages.  Rapid 
radial  expansion  occurs  at  early  times  in  the  central  region  r  <  1,  initially  at  the  inital  local 
sound  speed  c(r  —  0,  f  —  0),  but  later  slowing  down  to  the  ambient  sound  speed.  In  the 
second  stage,  expansion  occurs  at  a  velocity  only  moderately  higher  than  the  ambient  sound 
speed,  over  a  region  corresponding  to  the  final  channel  size,  i.e.,  a  low  Mach  number  shock 
propagates  over  this  region.  In  the  third  stage,  the  channel  returns  very  slowly  from  its 
minimum  density  state  to  the  final  equilibrium.  This  evolution  may  be  seen  in  Figs.  2. 1-2.8. 
The  data  is  also  replotted  in  compact  form  in  Figs.  2.14  and  2.15,  where  the  time  dependence 
of  p(r  —  0),  and  of  the  shock  radius  Rs(t)  are  shown. 
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We  note  from  Fig.  2.14  a  very  simple  rule  that  essentially  determines  the  end  of  the  first 
stage,  and  tells  us  how  long  it  takes  for  the  central  channel  to  be  created:  the  change  in  density 
is  90%  complete*  when*  t  =  0.7,  independent  of  P0.  It  takes  much  longer  for  the  remaining 
density  reduction  to  be  accomplished,  essentially  because  expansion  over  a  larger  volume  is  in¬ 
volved,  and  the  overpressure  becomes  weak. 

The  second  stage  can  be  said  to  terminate  when  minimum  density  is  reached;  this  time  fm 
increases  with  P0  because  the  size  of  the  low-density  channel  increases,  and  is  given  accurately 
(to  within  <  10%  for  all  values  of  P0  studied)  by  the  simple  formula  rm  -  1.2  x  ( P0  +  1)1/2. 

The  duration  of  the  third  stage,  i.e.,  the  return  to  pressure  equilibrium,  is  essentially  the 
time  for  sound  to  propagate  at  the  ambient  rate  over  the  reduced  pressure  region.  Typically  the 
return  to  equilibrium  has  an  e-folding  time  of  about  1.5/m. 

HI.  EXPANSION  FOLLOWING  ENERGY  DEPOSITION  IN  A 
PRE-EXCITING  CHANNEL 

A.  General  Remarks 

If  the  channel  is  heated  repeatedly  by  a  series  of  energy  pulses,  and  if  the  separation 
between  pulses  is  so  long  that  the  channel  essentially  comes  to  rest  between  pulses,  in  pressure 
balance  with  the  ambient  gas,  then  the  hydrodynamic  expansion  process  is  an  iterative  one,  in 
which  the  final  channel  state,  after  the  previous  pulse,  consitutes  the  initial  density  profile  for 
the  next  stage  of  expansion.  The  first  pulse  is  uniquely  simple  in  that  it  deposits  its  energy  in 
homogeneous  gas;  this  makes  possible  the  unambiguous  choice  of  dimensionless  coordinates  as 

•  But  note  that  if  ~P0  is  very  large,  then  the  final  10%  of  the  density  reduction  results  in  further  reduction  of  density  by 
a  significant  factor,  e.g.,  if  p(r  —  0.1,  the  factor  is  2. 
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in  Eqs.  (2.4,  2.5),  based  on  the  spatial  scale  a  and  the  time  scale  a/cA.  For  a  pulse  depositing 
energy  in  a  proformed  channel,  the  subsequent  hydrodynamic  expansion  depends  on  both  the 
beam  and  channel  radii,  as  well  as  on  the  overpressure;  furthermore,  the  time  scales  depend  on 
the  local  sound  speed  in  the  pre-existing  heated,  low-density  channel,  as  well  as  on  the  ambient 
sound  speed.  Thus  it  is  no  longer  possible  to  reduce  the  hydrodynamic  behavior  to  a  one- 
parameter  problem  by  choosing  an  obvious  set  of  dimensionless  coordinates.  Nevertheless,  we 
shall  see  that  the  calculation  of  the  channel  configuration,  after  it  comes  to  rest  in  pressure  bal¬ 
ance,  can  be  extended  to  the  present  case,  and  much  of  the  analysis  of  shock  properties,  pres¬ 
sure  undershoot,  and  time  scales  can  be  extended  in  a  qualitative  way.  We  also  consider  cases 
in  which  the  pulse  separation  is  such  that  the  channel  has  reached  a  uniform  pressure  which  is 
near  (but  not  quite  at)  ambient  pressure  when  the  next  pulse  arrives. 

We  continue  to  assume  for  specificity  that  the  overpressure  caused  by  the  n'h  pulse  in  a 
pre-existing  channel  has  a  Bennett  profile,  with  Bennett  radius  a.  As  in  Sec.  2,  Eqs.  (2.4),  we 
adopt  dimensionless  units  in  which  radial  distances  are  scaled  to  a,  times  to  a/cA  (where  cA  is 
the  ambient  sound  speed),  and  densities,  pressures,  and  temperatures  top,,,  PA,  and  TA.  Such 
dimensionless  quantities  are  identified  by  tildes.  In  addition,  it  is  now  convenient  to  define 
quantities  referred  specifically  to  the  time  of  arrival  of  the  n'h  pulse.  We  define  ~t„  =  t  -  f  „ 
as  the  dimensionless  time,  so  that  ~tn  =  0  at  the  time  of  arrival  of  the  n'h  pulse.  Also,  we  define 
r„_]  (  scaled  to  a)  as  the  dimensionless  coordinate  of  a  particular  element  fluid  at  t„  =  0.  As 
with  r0  in  Sec.  2,  r„_x  is  used  as  a  Lagrangian  coordinate,  i.e.,  a  label  for  an  element  of  gas.  In 
the  case  of  a  series  of  beam  pulses,  one  can  use  r0  as  the  Lagrangian  coordinate  while  following 
the  hydrodynamics  after  the  first  pulse,  r,  as  the  Lagrangian  coordinate  for  the  same  fluid  ele¬ 
ment  in  following  the  hydrodynamics  after  the  second  pulse,  etc.  The  Eulerian  coordinate  r, 
the  position  of  the  fluid  element  at  time  tn ,  can  be  regarded  as  a  function  of  and  By 
definition, 
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rn  =  r (?„_!,  ?«  ?„).  (3.1) 

We  also  define  Tn^l(r),  ?„_](?),  p„_i(r),  P„-\(r)  as  the  local  (dimensionless)  temperature, 

sound  speed,  density  and  pressure  characterizing  the  channel  at  t„  -  0~,  just  before  the  arrival 

of  the  n'h  pulse.  By  assumption, 

£„_](/■)- P„_i  (constant).  (3.2) 

If  the  pulse  spacing  is  long  enough  for  the  channel  to  come  to  rest  in  pressure  balance, 

P„  -  1  for  all  n,  (3.3a) 

and 


rn  =  r (rn_i,  ln  -  oo). 

For  the  special  case  of  the  first  pulse,  »  —  1  —  0, 


(3.3b) 


T„-i(r)  =  c„_,(r)  -  p„_i(r)  -  1. 

The  pressure  at  time  (f„  -  0+),  just  after  the  arrival  of  the  n'h  pulse,  is  taken  to  be 


(3.4) 


£(?„_„  t„  -  0+)  - 


1  + 


(1  +  V,)J 


(3.5) 


We  have  studied  a  wide  variety  of  cases  with  hydrodynamic  code  simulations.  Since  our 
purpose  here  is  to  elucidate  the  principal  features  of  the  hydrodynamic  evolution  and  arrive  at 
scaling  laws,  rather  than  to  overwhelm  the  reader  with  data,  we  show  only  the  results  of  a  few 
such  simulations  for  a  second  pulse  propagating  in  the  channel  formed  by  the  first  pulse,  and 
arriving  at  a  time  when  the  channel  is  near  its  minimum  density/pressure  state. 


B.  Pulse-to-Pulse  Evolution  of  Channel  Density 


We  consider  first  the  case  in  which  pulse  spacing  is  long  enough  so  that  the  channel 
comes  to  rest  in  pressure  balance  with  the  ambient  gas.  Following  exactly  the  same  reasoning 
that  led  to  Eqs.  (2.10c),  (2.12),  but  translating  it  into  the  notation  just  described,  we  find  that 
the  channel  profile,  when  it  comes  to  rest  after  n  beam  pulses,  can  be  expressed  in  terms  of  the 
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(3.7) 


We  recall  from  Sec.  2  that  exp  (— A5n(rn_|)/cp]  represents  a  small  (typically  <  20%) 
correction,  due  to  the  entropy  increase  AS„  (r„_|)  when  the  outward  propagating  shock  over¬ 
runs  the  fluid  element  defined  by  Lagrangian  coordinate  r„-\.  To  the  accuracy  required  to  cal¬ 
culate  channel  densities,  we  find  that  the  heuristic  formula  for  shock  strength,  Eq.  (2.23),  is 
also  applicable  to  the  present  case  of  a  pulse  depositing  energy  in  a  pre-existing  channel.  This 
is  evidently  true  because  the  shock  evolution  depends  principally  on  the  ratio  of  energy  deposit¬ 
ed  as  overpressure  to  energy  density  in  the  gas  before  the  passage  of  the  pulse;  the  latter  is  un¬ 
changed  by  the  presence  of  the  pre-existing  channel,  since  we  have  assumed  that  the  channel  is 
at  ambient  pressure.  Rewriting  Eqs.  (2.23,2.24),  we  have 


exp[-AS„/  cp\ 

where  the  density  jump  at  the  shock  is 


-  ,  +  5/7 

Pn  0  ~  Pn/Pn 

Pn  ^Pn/Pn  ~  1 


(3.8) 


(l,  Vi  <  0.8 

0„)  =  1  +  0.33  0.8  <3  (3.9) 

Pn  (1  +  1.7  3<V, 

Eqs.  (3. 6-3. 9)  now  give  a  complete  description  of  the  channel  density  profile,  when  pressure 
balance  is  restored  after  the  passage  of  a  pulse  in  a  channel,  in  terms  of  the  overpressure  Pm 
and  the  profile  of  the  pre-existing  channel.  By  iterating  these  equations  n  times,  we  obtain  the 
channel  profile  after  the  passage  of  n  pulses  in  initially  uniform  ambient  gas. 
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As  in  the  case  of  an  initial  pulse  depositing  its  energy  in  uniform  gas,  we  find  that  for  ail 
cases  with  central  overpressure  Pon  >  1,  the  pressure  profile  P(r,  7„)  reaches  a  very  broad, 
shallow  minimum  value  Pm,  slightly  less  than  unity,  at  sometime  Jm„.  This  is  shown  for  several 
fluid  code  runs  in  Figs.  3.4-3.6.  In  fact,  the  value  of  Pm  is  found  to  depend  only  on  Pon,  and 
not  on  the  channel  properties,  i.e.,  it  is  given  by  Fig.  2.13,  even  in  the  presence  of  a  pre¬ 
existing  channel.  For  practical  purposes,  it  is  sufficient  to  estimate  P„  -  0.88  for  all  cases  with 
Kn  >  I- 


It  might  be  desirable,  in  some  cases,  to  choose  the  pulse  spacing  such  that  each  fresh 
beam  pulse  arrives  at  the  time  ~tm„  when  the  channel  is  deepest,  rather  than  after  it  has  returned 
to  pressure  balance.  In  this  case,  since  P{r)  =  Pm  ~  0.88  is  essentially  constant  over  the 
channel  (and  the  same  for  each  successive  pulse),  Eqs.  (3.6-3. 9)  remain  formally  unchanged 
for  n  >  2.  We  note,  however,  that  the  overpressure  parameter  Pon  which  occurs  in  Eq.  (3.6), 
is  defined  in  Eq.  (3.5)  relative  to  the  channel  pressure  P„-\  at  t„  -  0;  thus,  for  a  given  amount 
of  energy  deposition  by  the  n'h  pulse,  Pon  is  larger  if  the  channel  is  at  minimum  pressure  rather 
than  ambient  pressure.  The  evolution  following  the  first  pulse  begins  with  P{r)  -  1  and  ends 
with  P(r)  ~  0.88,  so  for  the  first  pulse  Eq.  (3.6)  is  modified  to 


Pi  (?.) 


1  + 


r0i 


(i  +  f0y 


-i  h 


(0.88) exp  (-AS(?b_,)/c/>], 


(3.10) 


as  discussed  in  detail  in  Sec.  2C. 


In  practice,  the  dependence  of  the  channel  evolution  on  the  exact  pulse  spacing  is  rather 
weak,  and  interpolation  is  easily  performed,  provided  the  channel  pressure  is  anywhere  close  to 
the  ambient  pressure  when  the  next  pulse  arrives. 
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C.  Time  Evolution 

The  time  scales  for  radial  expansion,  following  energy  deposition  by  a  beam  pulse  in  a 
preformed  channel,  are  analogous  to  those  for  energy  deposition  in  uniform  ambient  gas,  dis¬ 
cussed  in  Sec.  2D.  The  situation  is  more  complicated  because  the  sound  speed  c„_i(r)  in  the 
pre-existing  channel  is  now  a  function  of  r. 

The  time  evolution  of  p(r  —  0 ,tn)  and  of  the  shock  radius  ^s(r„)  are  shown  in  Figs.  3.10 
and  3.11,  for  three  cases  run  on  the  fluid  code.  Once  again,  we  note,  in  three  figures  and  in 
Figs.  3.1-6,  the  existence  of  three  stages  of  development.  In  the  first  stage,  rapid  expansion 
occurs  over  r  <  1,  initially  at  the  sound  speed  c(r=  0,  —  0+)  of  the  heated  central  region. 

Since  the  pre-existing  channel  is  typically  much  broader  then  r  —  1,  this  stage  samples  only  the 
uniform  central  region  of  the  channel.  Thus  this  stage  of  evolution  is  similar  to  that  which 
occurs  for  energy  deposition  in  uniform  ambient  gas,  except  that  the  ambient  sound  speed  is 
replaced  by  the  sound  speed  c„ _ j ( r )  in  the  pre-existing  channel.  Indeed  we  see  in  Fig.  3.10 
that  the  change  in  density  at  f  —  0  is  90%  completed  when  t„  =  0.7  in  analogy  with  the 

conclusion  of  Sec.  2D  for  the  first  pulse. 

In  the  second  stage,  a  low  Mach  number  shock  propagates  over  the  final  channel  radius, 
at  a  speed  moderately  higher  than  the  local  sound  speed  c„_  j  (?)  in  the  pre-existing  channel. 
Since  c„_ i  (? )  varies  widely  over  the  channel,  the  time  scale  for  this  process  depends  on  both 
the  strength  of  the  n'h  pulse  and  the  depth  and  width  of  the  pre-existing  channel,  making  it 
difficult  to  give  any  simple,  general  formula  for  this  time  scale.  In  general,  a  minimum  pres¬ 
sure  is  reached,  which  is  nearly  constant  over  a  broad  region,  and  this  occurs  at  a  time  tm„ 
which  varies  from 

L  =  (1  to  2)  x  ( Pon  +  1),/2 
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for  a  pulse  which  is  strong  enough  to  significantly  broaden  the  channel  to 

In  ~  1.2  X  (Pon  +  l)'/2  f-J{2  (r  -  0) 

for  a  pulse  which  is  too  weak  to  produce  much  channel  broadening.  However,  the  minimum 
pressure  itself  depends  essentially  only  on  P0„,  and  not  on  the  channel  properties,  and  is  given 
to  good  accuracy  by  Fig.  2.13;  for  the  whole  range  Pgn  ^  1,  this  minimum  pressure  is 

Pnm  -  0.85  to  0.9. 

In  the  third  stage,  the  channel  slowly  returns  to  pressure  equilibrium  at  P{r)  =  1.  As  in 
the  case  of  a  first  pulse  into  ambient  gas,  this  is  essentially  a  linear  sound  wave  propagation  pro¬ 
cess,  and  the  e-folding  time  for  the  approach  to  equilibrium  is  typically  about  1.5/m„. 

IV.  WEAK  BEAM  PULSE  -  LINEARIZED  HYDRODYNAMICS 

In  this  report  we  have  studied  the  nonlinear  hydrodynamic  response  to  a  beam  pulse 
which  instantaneously  heats  the  gas  it  traverses.  If  the  pulse  happens  to  be  weak,  i.e.,  the  over¬ 
pressure  is  much  smaller  than  the  ambient  pressure  PA,  and  the  air  is  initially  uniform,  then 
the  hydrodynamics  can  be  linearized  and  solved  analytically  Such  calculations  are  of  direct 
interest  in  application  to  a  weak  pulse  or  a  series  of  weak  pulses,  and  in  addition  provide  valu¬ 
able  insight  into  the  hydrodynamic  features  of  the  general  problem.  Results  of  this  type  have 
been  presented  by  Fader"  for  the  special  case  of  a  Gaussian  radial  heating  profile.  We  extend 
these  results  somewhat  here,  and  treat  in  addition  the  cases  of  Bennett  and  step  function  heat¬ 
ing  profiles. 

In  the  notation  of  this  section,  ambient  quantities,  e.g.  pressure  PA,  density  pA,  sound 
speed  cA ,  temperature  TA,  are  designated  by  subscript  A.  Perturbed  quantities  are  defined 
A P(r,t),  &p(r,t),  \c(r,t),  &T(r,t),  and  fluid  radial  velocity  A u(r)  and  are  all  assumed  to  be 
small,  of  first  order.  The  problem  becomes  essentially  one  of  linear  sound  propagation  driven 
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by  a  given  initial  overpressure.  Weak  shock  waves  do,  strictly  speaking,  eventually  form,  but 
this  occurs  at  a  very  late  time,  and  any  accompanying  entropy  increases  are  higher  order,  and 
thus  negligible.  Thus  the  isentropic  solution  of  Eqs.  (2.10c),  (2.12)  for  the  density  profile 
holds  at  late  times  (when  pressure  balance  is  restored).  Furthermore,  these  equations  reduce 
to  the  very  simple,  linearized  form 


A p(r,t  =  °o)  =  -  pA 


\P(r,0+) 


&P(r,  0+) 


A  T(r,t  =  oo)  = 


A  T(r,  0+) 


The  late-time  density  and  temperature  profiles  have  exactly  the  shape  and  width  of  the  in¬ 
itial  overpressure  profile.  We  see  from  Eq.  (4.2)  that  a  fraction  y-1  of  the  inital  overpressure 
thermal  energy  remains  in  the  channel  when  it  comes  to  rest.  The  remainder  of  this  energy  is 
carried  to  r  =  °°  by  an  outward  propagating  acoustic  pulse. 

The  time  dependent  hydrodynamic  process  is  determined  by  the  linearized  fluid  equations, 

^ +  t "h  <,a")  -  °-  <«> 

with  the  linearized  adiabatic  equation  of  state 

dP  HP(r,t)  -  Af(r.O)  2  , 

dp  &p(r,t)  °A’  (4.5a) 


A  Hr./)-  ±I&2l+  1-1  Mkdl 


and  the  initial  conditions 


Ap  =  0,  \P(r,  0)  given. 
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We  note  that  the  /c-integration  in  (4.15b)  can  be  performed  in  closed  form,  thus  reducing 
the  solution  to  a  single  quadrature.  However  we  take  a  different  route.  For  computational  pur¬ 
poses,  it  is  convenient  to  use  the  relation 

J0(kr)  -  —  f  dx  sin(/cx)  Or2  -  r2)~V2  (4.16) 

n  Jr 

to  recast  Eq.  (4.15a)  in  the  form 

t)  «  —  (irc^p^)  §  dx(x2  -  r2)~l/2 

x  J*o  dk  {cos(/c(cr  -  x)]  -  cos(/c(cr  +  jc)]}  APk(ff).  (4.17) 

The  /c-integration  in  (4.17)  can  be  peformed  analytically  for  many  choices  of  AP(r,  0). 


Next  we  consider  three  particular  choices  of  AP(r,  0):  the  Bennett  profiie, 

,  AP0 

AP{r0)-—--2-/a2)2, 

for  which  (with  K\  the  modified  Hankel  function), 

APk( 0)  =  |  a2k  AP0  K\(ka)\ 

the  step  function. 


(4.18a) 


(4.18b) 


for  which 


AP(r,  0)  -  AP0  H(a  -  r)  = 


A  P,  r  ^  a 
0,  r  >  a 


APk( 0)-  AP0  J\(ka)\ 

and  the  Gaussian  profiie, 

AP(r,  0)  =  AP0  exp(-r2/a2), 

for  which 


(4.19a) 


(4.19b) 


(4.20a) 


APk(0)  -  j-a2  exp 


~7  ><2a2 

4 


(4.20b) 


24 


NRL  MEMORANDUM  REPORT  4073 


In  all  three  cases,  the  /e-integration  in  Eq.  (4.17)  can  be  performed  explicitly,  giving  the  follow¬ 
ing  results  for  0  in  terms  of  single  quadratures  (using  a  further  transformation,  x  —  r  cosh  z  to 
remove  the  logarithmic  singularity  at  r  -  x).  For  a  Bennett  profile, 

f)  ~  ~^PaC~a  X  dz  (l^'  “  r  coshz)2  +  fl2J  3  2 

-[(c^r  +  r  coshz)2  +  a2J  V2|;  (4.21) 

for  a  step  function  profile 

cAt  -  r  coshz 1 2|  1 


t)  -  - 


A  P0a 


irpAcA  Jo 


/; 


dz 


x  H(a  -  | cAt  -  r  coshzl) 


cAt  +  r  coshz 

1 

2 

1  - 

a 

H(a  -  Ic^t  +  r  coshz |) 

and  for  Gaussian  profile. 


0(r,  r)  - 


-A  P0a 
2 PaCa^'2 


f"dz 

|exp 

_ 

ct  -  r  coshz 

2 

J  0 

a 

-exp 


ct  +  r  coshz 
a 


(4.22) 


(4.23) 


The  single  integrals  in  Eqs.  (4.21)-(4.23)  have  been  carried  out  numerically.  The  quanti¬ 
ties  \P(r,t),  bp(r,  t),  A T(r,t)  are  determined  from  0(r,/)  via  Eqs.  (4.5).  The  results,  in 
each  case,  are  universal  when  expressed  in  ten.-*  of  the  dimensionless  variables 
z  “  r!a’  1  “  cAtl°i  P  =  A P/PA,  p'«*  A p/pA,  7"'=  \T/Ta.  In  Figs.  4.1  through  4.5  we  show 
profiles  of  A p(r)  and  A P(r)  at  a  sequence  of  times,  for  the  Bennett  and  step  function  profiles, 
and  in  Figs.  4.5  and  4.6  the  time  dependence  of  p(r  -  0)  is  displayed. 
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In  each  case  an  outward  propagating  sound  pulse,  and  an  inward  propagating  rarefacton 
wave,  originate  at  the  point  of  maximal  Id/Mr,  o)/3r|,  both  propagating  at  cA.  This  is  particu¬ 
larly  clear  for  the  step  function  initial  profile,  for  which  there  is  no  disturbance  at  a  point  r  until 
t  >  \a  -  r  \  / cA.  The  outward  pulse  propagates  to  r  =  °°,  carrying  with  it  a  constant  fraction 
(1  -  y-1)  of  the  energy  deposited  in  the  medium.  The  density  and  pressure  at  points  near  the 
axis  drops  monotonically  until  the  time  ~tm  when  the  rarefaction  pulse  hits  the  axis.  At  this 
time,  the  pressure  on  axis  reaches  a  minimum  Pm  ( r  —  0)  and  then  rebounds  slowly  and  mono¬ 
tonically  until  it  reaches  PA.  There  is  no  perceptible  series  of  oscillations  or  "bounces"  of  the 
channel  that  is  left  behind.  The  value  of  the  density  minimum,  pm(r  =  0)  =  p(r  =  0,  t  —  7m), 
is  given  in  Table  4  for  each  of  the  three  assumed  profiles  of  P{r0),  step  function,  Gaussian, 
and  Bennett. 

All  of  these  features  carry  over,  with  obvious  modification,  to  the  nonlinear  hydrodynam¬ 
ics  studies  in  Sec  2.  However  one  difference  is  that  in  the  nonlinear  case,  with 
i\P(0,  0)  >  PA ,  we  found  that  at  time  tnm  the  pressure  is  practically  constant,  at  its  minimum 
value,  over  the  reduced-density  channel.  In  the  present  linear  regime,  the  pressure  wave  of 
course  has  the  same  characteristic  wavelength  as  the  density  wave;  thus  the  pressure  profile  at 
tm  is  no  broader  than  the  density  channel,  and  different  parts  of  the  channel  reach  minimum 
pressure  at  different  times. 

V.  CONCLUSIONS 

In  this  report,  we  have  studied  the  hydrodynamics  of  channel  formation  over  a  wide  range 
of  parameters  ranging  from  very  weak  overpressures  up  to  overpressures  ~40  times  ambient. 
Channel  cooling  and  variation  of  the  adiabatic  index  y  with  p  and  T  have  been  omitted.  The 
main  conclusions,  for  the  case  of  heating  pulses  with  smooth  (Bennett)  radial  profiles,  are  as 
follows. 
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(i)  The  density  reduction  at  the  center  of  the  channel,  when  the  channel  motion  comes 
to  rest,  is  easily  and  accurately  calculated  by  the  isentropic  model,  for  any  number  of  widely 
spaced  pulses.  Results  are  given  in  Table  3  and  Eq.  (2.12). 

(ii)  The  isentropic  model  of  Table  3  and  Eqs.  (2.10c),  (2.12)  is  accurate  to  within 
<  20%  per  pulse  in  calculating  the  overall  shape  and  width  of  the  channel,  for  overpressure 
P0  <  40  (it  is  more  accurate  for  smaller  overpressures). 

(iii)  An  adiabatic  model,  extended  to  include  shock  wave  heating,  is  highly  accurate  in 
calculating  overall  channel  shape  and  width.  Results  are  given  in  Table  3  and  Eqs.  (2.8,  10c, 
23,  24),  (3.6-9). 

(iv)  During  the  hydrodynamic  evolution,  the  channel  undershoots  its  final  density  by  a 
factor  — •  1 0% ,  essentially  constant  over  the  channel,  and  practically  independent  of  P0  for 
P0  >  1  It  then  returns  slowly  to  pressure  balance,  without  "ringing". 

(v)  The  time  scales  for  channel  formation  near  r  =  0,  for  the  whole  channel  to  reach 
minimum  density,  and  for  the  return  to  pressure  equilibrium,  are  distinct,  and  are  given  in 
Table  2.  The  first  time  scale  is  typically  up  to  an  order  of  magnitude  shorter  than  the  other 
two. 

(vi)  For  very  weak  beam  pulses,  the  complete  time-dependent  linearized  hydrodynamics 
can  be  solved,  in  terms  of  quadratures.  Results  are  given  in  Sec.  4. 

(vii)  The  self-similar  blast  wave  theories  of  classical  hydrodynamic  lore5  are  only  applica¬ 
ble  to  beam  heating  much  stronger  than  are  considered  here. 
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Table  I  —  Notation 


a 

Typical  (Bennettl  radius  of  beam  and  of 
beam-heated  region 

cn 

Specific  heat  at  constant  pressure 

c 

Specific  heat  at  constant  volume 

y 

=  Lp/c, ,  adiabatic  index,  taken  to  be  — 

P(r.  1 ) 

Pressure 

p(r,  1 ) 

Density 

T(r,  i ) 

Temperature 

c(r,  l ) 

Sound  speed 

P, 

Ambient  pressure 

Pi 

Ambient  density 

T, 

Ambient  temperature 

Ci 

Ambient  sound  speed 

r 

Position  coordinate  (Fulerian) 

1 

Time,  measured  since  arrival  of  first  beam  pulse 

r„ 

Time,  measured  since  the  nth  beam  pulse  (r,  =  t) 

T„ 

Time  of  arrival  of  nth  pulse,  r,  =  0. 

r„ 

Initial  position  of  a  fluid  element,  used  as  a 

Lagrangian  coordinate,  i.e.  a  label  for  that  element 

rn  -i 

Position  of  a  fluid  element  at  t  —  r„, 
arrival  time  of  nth  pulse.  Also  used  as  a 

Lagrangian  coordinate 

P„i(r) 

Pressure  at  time  t  =  r„ 

Pn-\lr) 

Density  at  time  t  -  t„ 

T.,  ,(/•) 

Temperature  at  time  i  =  r„ 

C„  ,(r) 

Sound  speed  at  time  i  =  r„ 

P,Jr) 

Overpressue  due  to  energy  deposition  by  nth  pulse 

P.,„ 

P„  (r-o)/P,l  ,,  used  to  label  strength  of  nth  pulse 

r,„ 

Lagrangian  coordinate  of  fluid  element  where  shock  forms  1 

A 

Time  when  shock  forms 

«,(/» 

Shock  location  at  time  i 

c, 

Time  when  channel  reaches  minimum  pressure 

A.V(r„ ... 

Lntropy  increase,  due  to  shock 

A  =  P/P , 

Dimensionless  forms  are  defined  similarly 

P  =  P/Pt 

for  various  subscripted  forms 

f  =  r/r, 

c  =  c/c., 
r  =  r/a 
r  =  Cfl/a 

of  P,  p.  7,  c,  /  and  / 
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Table  2  —  Time  Scale  for  Hydrodynamics 
1)  Density  reduction  at  r  =  0. 


p(r  =  0)  reduced  90%  of  the  way 
to  its  asymptotic  value,  i.e. 

p(0,  7„)  =  p(0,  °o)  +0.1  (p„_,  (0)  -  p( 0,  °°)] 

■ _ _ _ _ _ 

2)  Minimum  density  and  pressure,  entire  channel  r„— 1„ 

First  pulse 

Strong  pulse  in  existing  channel 
Weak  pulse  in  existing  channel 

3)  Channel  reaches  pressure  balance  with  ambient  gas 


l  -0.7  If,,.,  (0)]-"2 


»«■  -  12(1  +  P,,)'12 

r„„,  ~  (1  to  2)  (1  +  P„)'/2 

l„„  =  1.2  (1  +  P0)'nlt„ (0))' 

/«  3 


Table  3  —  Channel  Density  After  Gas  Reaches  Pressure  Balance 

(1)  Channel  density  at  r  =  0 

p„(0)  ”  p„_,  (0)  (1  +  Pnn)  _,/>  (lsentropic  expansion) 

(2)  Channel  density  profile 

We  label  an  element  of  gas  by  its  position  r„_,  before  the  pulse  arrives.  r„  is  the  position  of 
the  same  gas  element  after  pressure  balance  is  re-established. 


-2  C  "- 1  ji~  Pn-\  (/»-!  ) 

“ l,  d(r-')  )  ’ 


p..  ,  ,  P„„ 

T-  T  =  +  7,“  Z~2  o  exp[— AS„  (r„_,)/cp], 

p„_,  (/■„-!>  (1  +  ',,-1  )2 


-  I  /  _  +  5/7 

expl-AS„  ( _ j ) / c,, ]  =  ^  -p  p  (entropy 

P  6p  /p+  -  1 


jump  due  to  shock). 


1,  /*„  <  0.8 

—  (r„-|)  -  1  +  0.33  P„'„/4  /*„,  0.8  <  r  <  3(density  jump  at  shock). 

1  +  1.7  'V'/2.  3  < 


Final  coordinate  ?,  =  rt/oof  a  fluid  element  whose  initial  coordinate  is 
as  a  function  of  initial  central  overpressure  Pn. 
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Fig.  1.3  —  A  graphical  iterative  method  of  going  from  the  initial  (dimensionless)  coordinate  r0  of  a  fluid  ele¬ 
ment,  to  the  coordinate  >|  when  the  channel  has  returned  to  pressure  equilibrium  after  the  first  pulse,  to  the 
value  r}  at  pressure  equilibrium  after  the  second  pulse,  etc.  For  the  case  shown,  P0,  -  4.32,  fl02  “  34.9,  and 
r0-  1. 
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Fig.  2.1  —  Radial  profiles  of  (he  scaled  density  (p  =  p/p4)  disturbance,  for  a  Bennett  ini¬ 
tial  pressure  distribution  (of  radius  a)  with  initial  overpressure  P0  =  /y  P4  -  4.20  at  >  = 
r/a  -  09.  Subscript  A  denotes  ambient  values.  The  times  t  =  icja  -  0.33,  0.66,  and 
9.66,  respectively,  are  close  to  the  shock  onset  time,  the  time  when  the  shock  is  strongest 
and  a  time  when  the  channel  density  is  close  to  but  slightly  below  its  final  equilibrium 
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Fig.  2.2  —  Radial  profiles  of  the  scaled  density  p,  for  a  Bennett  initial  pressure  distribution  with 
initial  overpressure  ~Po  -  10.5  at  >  -  0.  Times  are  chosen  by  analogy  to  those  in  Fig.  2.1. 
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Radial  profiles  of  scaled  density  p,  for  a  Bennett  initial  pressure  distribution 
with  initial  overpressure  /n  —  0.1  at  r  •»  0. 
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Fig.  2.5  -  Radial  distribution  of  the  scaled  pressure  P  =  P/ P„  for  a  Bennett  initial  pressure  distribution  with  initial 
overpressure  ~P0  —  42.0  at  r  —  0.  The  times  f  —  0.33,  0.66  and  9.66  respectively  are  close  to  the  shock  onset  time,  the 
time  when  the  shock  is  strongest  and  a  time  when  the  channel  density  is  close  to  but  slightly  below  its  final  equilibrium 


Radial  distributions  of  the  scaled  pressure  P,  for  a  Bennett  initial  pressure  distribution 
with  initial  overpressure  P0  -  0.1  at  r  -  0.  Times  are  the  same  as  Fig.  2.4. 


Fig.  2.9  —  The  radial  density  profile  for  the  case  P—  42.0  as  given  by  the  isentropic  approximate  model  of  Eqs. 
(2.10c),  (2.12),  is  plotted  as  the  solid  curve.  The  crosses  show  the  density  profile,  as  given  by  the  complete  analytic 
model  of  Eqs.  (2.33),  (2.8),  and  (2.10c),  which  includes  shock  heating.  The  dashed  curve  shows  the  fluid  code  result 
for  the  same  case  at  late  time. 
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Fig  2.10  —  Similar  to  Fig.  2.9,  but  for  P0 


10.5. 
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Fig.  2.11  -  Similar  to  Fig.  2.9,  but  for  P0-  1.8.  Shock  heating  plays 
no  significant  role  in  this  case  of  weak  expansion. 
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Fig.  2.12  —  Similar  to  Fig.  2.11,  but  for  P0  -  0.1. 


45 


Fig  2.13  -  The  solid  curve  is  the  pressure  minimum  P„  s  P„/Pt  plotted  as  a  function  of  the  initial  overpressure 
A>  =  (/’(r  -  0,  /  -  0)//%)  -  1,  as  gi^en  by  fluid  code  results  for  a  first  pulse  in  uniform  ambient  gas.  The  crosses  are 
also  fluid  code  values  of  P„,  for  three  o.fTerent  cases  of  energy  deposition  by  a  pulse  in  a  preformed  channel. 
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Fig.  3.1  —  Radial  profiles  of  ihe  scaled  density  disturbance  (p  =  p/p4).  The  solid  curve  shows  the  density  profile  at 
time  1 1  -  9.66  following  the  passage  of  a  first  pulse  with  P0I  -  43.2.  The  other  curves  show  the  disturbance  at  various 
limes  i2  after  the  passage  of  a  second  pulse  with  overpressure  P02  -  34.9.  The  times  J2  ~  0.33,  40,  and  9.66  are  close 
to,  respectively,  the  shock  onset  time,  Ihe  time  when  the  central  pressure  reaches  the  ambient  value,  and  a  time  when 
the  channel  density  is  close  to  but  slightly  below  its  final  equilibrium  value 


4  -  Radial  profile  of  the  scaled  pressure  (P  =  P/Pt),  al  a  lime  just  before  the  passage  of  the  second 
pulse  (solid  curve),  and  at  various  times  after.  Taken  from  the  same  fluid  code  simulation  as  Fig.  3.1 
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Fig.  3.7  —  The  radial  profiles  of  the  scaled  density  p  due  to  the  first  and  the  second  pulses  are  shown  for  late  time  in 
dotted  curves.  The  solid  curves  show  the  iteration  results  for  two  pulses  obtained  by  the  isentropic  approximation. 
The  crosses  show  the  interation  results  for  two  pulses  obtained  from  the  complete  analytic  model,  including  the  shock 
entropy  production  given  by  Eq.  (3.9).  The  central  overpressures  produced  initially  by  the  two  pulses  are  P0|  -  43.2 
and  hi  -  34.9  (same  fluid  code  simulation  as  Figs.  3.1  and  3.4). 
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Fig  3  9  —  The  dolled  curves  show  ihe  fluid  code  density  profiles  al  late  limes,  due  lo  the  first  and  second  pulses  pro¬ 
ducing  overpressure  Pol  —  1  81  and  P02  —  1.6  The  solid  curves  show  the  iteration  results  obtained  by  the  isentropic 
approximation.  Shock  heating  plays  no  significant  role  in  these  weak  expansions. 
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Fig.  3.10  —  The  lime  evolution  of  the  central  density,  during  expansion  following  the  passage  of  a  second  pulse  in  a 
preformed  channel.  The  three  cases  of  Figs.  3.1,  3.2,  and  3.3  are  shown,  i.e.,  respectively  P02  —  34.9  (solid  curve),  8.9 
(long  dashed)  and  1.6  (short  dashed).  The  vertical  scale  is  chosen  so  that  the  initial  value  is  1.0  and  the  final  value  is 
zero  in  all  cases.  The  time  is  scaled  to  a/c  (0,0),  where  c(0.0)  is  the  sound  speed  at  r  -  0  just  before  the  second  pulse 
arrives.  The  zero  of  time  is  when  the  second  pulse  arrives. 
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Fig.  4.2  —  Radial  profiles  of  the  pressure  perturbation  it  a  sequence  of  times  (Bennett  initial  overpressure  profile). 
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Fig  4.4  -  Radial  profiles  of  the  pressure  perturbation  at  a  sequence  of  times  (step  function  initial  overpressure  profile) 
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Time  dependence  of  density  perturbation  and  pressure  perturbation  at  the  origin,  for  Bennett  initial 
overpressure  profile  (same  curve,  different  vertical  coordinates). 


